function res = predictiveRegression(Data,m,exoRegressor)

% The egression reads:
% rx_t+m,k = alpha + beta*(y_t_k-y_t_m) + eps_t+1
% where 
% rx_t+m,k = -(k-m)/4*y_t+m,k-m + k/4*y_t_k - 1/4*y_t_m
%
% INPUT: Must be a data matrix of dim T * numYields, 
% where yields have maturities of 1 quarter, 2 quarters,3 quarters, ...
[T,n]      = size(Data);
matSelect  = m:m:n;
shorRate   = Data(1:T-m,m);
yieldsUsed = Data(:,matSelect);
numYields  = size(yieldsUsed,2);
if nargin > 2 
    Betta      = NaN(2+size(exoRegressor,2),numYields);
    Betta_tstat= NaN(2+size(exoRegressor,2),numYields);
    Betta_se   = NaN(2+size(exoRegressor,2),numYields);
else
    Betta      = NaN(2,numYields);
    Betta_tstat= NaN(2,numYields);
    Betta_se   = NaN(2,numYields);
end
Ydata      = NaN(T-m,numYields);
Xdata      = NaN(T-m,numYields);
R2         = NaN(1,numYields);
for i=1:numYields-1
   % Note that m*i gives k-m in relation to the notation above
   Y  = (-m*i/4*yieldsUsed(1+m:T,i) + m*(i+1)/4*yieldsUsed(1:T-m,i+1) - m/4*shorRate);
   x  = (yieldsUsed(1:T-m,i+1)-shorRate);
   X  = [ones(T-m,1) x];
   if nargin > 2 
       X = [X exoRegressor(1:T-m,:)];
   end
   resOLS            = nwest(Y,X,m);
   Betta(:,i+1)      = resOLS.beta;
   Betta_tstat(:,i+1)= resOLS.tstat;
   Betta_se(:,i+1)   = resOLS.beta./resOLS.tstat;
   Ydata(:,i+1)      = Y;
   Xdata(:,i+1)      = x;
   R2(1,i+1)         = resOLS.rsqr;
end
% Output
res.betta       = Betta;
res.betta_tstat = Betta_tstat;
res.maturities  = matSelect;
res.betta_se    = Betta_se;
res.Ydata       = Ydata;
res.Xdata       = Xdata;
res.R2          = R2;

